function [smooth_st, smooth_counter, countmat, counter_countmat] = lmj_counterfactual(funcmod,  par1, par2, Y, trainvec, varargin)
% (1) Call LMJ_KFITERCENTRAL using par1
[stt,etamat,smooth_st,yfor,vstar,logL,countmat,countobs,countszero]=lmj_kfiltercentral(par1,funcmod,Y,trainvec, varargin{:}); %,solveopt,addsol); 
[T, ns, nx] = size(countmat); %number of shocks in total; 

% check smooth_st and countmat is correctly generated
check1 = nan(nx, 1); 
for ii = 1:nx; 
    check1(ii) = max(abs(sum(countmat(:, ii, :), 3) - smooth_st(:, ii))); 
end
if ~all(check1<1e-8); 
    error('variable countmat and/or smooth_st not defined corrected')
end

% (2) solve model using par2
[GG2, RR2, CONS2, eu2, SDX2, ZZ2]=feval(funcmod,par2,varargin{:});
smooth_counter = nan(T, ns); 


smooth_counter(1, :) = sum(countszero, 2); % initialization
for tt = 2:T; 
    smooth_counter(tt, :) = (GG2*smooth_counter(tt-1, :)'+ RR2*etamat(tt, :)')'; 
end

counter_countmat = nan(T, ns, nx); 
for ii = 1:nx; 
    tempeta = zeros(nx, T);
    tempeta(ii, :) = etamat(:, ii)';
    stsim = nan(ns, T); 
    stsim(:, 1) = countszero(:, ii);
   
   
    for tt = 2:T; 
        stsim(:, tt) = GG2*stsim(:, tt-1) + RR2*tempeta(:, tt); 
    end
    counter_countmat(:, :, ii) = stsim(:, :)';
end

% check smooth_counter and counter_countmat is correctly generated
check2 = nan(nx, 1); 
for ii = 1:nx; 
    check2(ii) = max(abs(sum(counter_countmat(:, ii, :), 3) - smooth_counter(:, ii))); 
end    
if ~all(check2<1e-8); 
    error('variable countmat and/or smooth_st not defined corrected')
end


end 